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Abstract. We present a new global fit of the 19-dimensional phenomenological Minimal 
Supersymmetric Standard Model (pMSSM-19) that comply with all the latest experimental 
results from dark matter indirect, direct and accelerator dark matter searches. We show that 
the model provides a satisfactory explanation of the excess of gamma-rays from the Galactic 
centre observed by the Fermi Large Area Telescope, assuming that it is produced by the 
annihilation of neutralinos in the Milky Way halo. We identify two regions that pass all the 
constraints: the first corresponds to neutralinos with a mass ~ 80—100 GeV annihilating into 
WW with a branching ratio of 95% ; the second to heavier neutralinos, with mass ~ 180 — 200 
GeV annihilating into tt with a branching ratio of 87%. We show that neutralinos compatible 
with the Galactic centre GeV excess will soon be within the reach of LHC run-II - notably 
through searches for charginos and neutralinos, squarks and light srnuons - and of XenonlT, 
thanks to its unprecedented sensitivity to spin-dependent cross-section off neutrons. 
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1 Introduction 

The existence of dark matter (DM) is now robustly established [1-4] and its cosmologi¬ 
cal abundance measured with high precision [5]. Yet, the fundamental nature of the most 
abundant matter component in the Universe is unknown. According to the most promising 
theories, DM is a new fundamental particle. As a consequence, the search for DM is also 
a search for new physics beyond that of the well-known elementary particles, as laid down 
in the Standard Model (SM). Weakly interacting massive particles (WIMPs) are the leading 
DM candidates: they arise in many extensions of the Standard Model of particle physics, and 
naturally achieve the appropriate relic density through self-annihilation the early Universe. 
WIMPs can be searched for with three detection strategies: direct detection of the energy 
recoil of nuclei scattering off DM particles; indirect detection of the final stable products 
of DM annihilation or decay, as for example gamma rays; and accelerator searches for new 
particles, in particular at the Large Hadron Collider (LHC). 

An excess in the 7 -ray emission from the centre of our Galaxy has been discovered in 
data from the Large Area telescope (LAT), aboard the Fermi satellite [6-16]. The nature 
of the so-called Fermi GeV excess remains a mystery. Several explanations have been put 
forward, the most exciting of which is perhaps DM annihilation in the halo of the Milky Way 
(see for example [13, 14, 17]). Astrophysical processes have also been suggested: the emission 
from a population of dim unresolved sources [18-21], and the inverse Compton emission from 
a new population of cosnric-ray, either from time-dependent events taking place in the past 
of the Galaxy [22-24] or from the high star formation activity in the inner Galaxy [25]. 

Very recently, two reanalyses of the gamma-ray emission from the inner Galaxy found 
strong evidence for the Fermi GeV excess being due to hundreds to thousands of dim or 
unresolved point sources [26, 27], most likely millisecond pulsars [18]. Even more recently, it 
was shown that part or all of the required millisecond pulsar population could originate from 
the disruption of globular clusters by tidal forces in the inner Galaxy [28]. Directly detecting 
some of the millisecond pulsars in the inner Galaxy by radio observations is the next critical 
step for fully establishing this scenario [26]. 
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It is however striking that the Fermi GeV excess spectrum and spatial distribution are 
well fitted by what is expected from DM annihilation. The excess could be the first non- 
gravitational signal of DM particles. It is thus urgent to either corroborate or disprove the 
DM particle nature of the Fermi GeV excess in the framework of concrete models for physics 
beyond the SM. Supersymmetry (SUSY) is one of the best-motivated classes of renormal- 
izable extensions of the SM, which can accommodate a stable DM particle together with 
new degrees of freedom that mediate interactions. The most generic R-parity conserving 
and phenomenologically viable supersymmetric model is the phenomenological Minimal Su¬ 
persymmetric Standard Model (pMSSM) [29]. In this work we address the following key 
question: can models of the pMSSM explain the observed properties of the Fermi GeV ex¬ 
cess, while retaining consistency with other experimental data? And if so, what are the 
detection prospects for future direct detection and collider experiments? 

Until now, MSSM scenarios [30] could not reproduce the Fermi GeV excess as observed 
by [14], as it was impossible to obtain in this framework a light neutralino (m x ~ 30 - 40 
GeV), as required to fit the Fermi GeV excess spectrum, which could also account for the 
cosmological DM as measured by Planck [5]. However, [15] demonstrated that higher WIMP 
masses and annihilation channels different from b-quark pairs can give a good fit to the 
Fermi excess, owing to the freedom allowed by background model systematics. By fitting the 
GeV excess data of [15], it has been shown that viable solutions in the MSSM exist [31]. In 
the context of the pMSSM, [32] demonstrated that a re-assessment of the theoretical uncer¬ 
tainties in the DM signal spectra opens up a new phenomenology at the LHC experiments. 
Here we present the first systematic study of the pMSSM parameter space through global 
fitting techniques. This approach exhaustively covers all possible phenomenological signa¬ 
tures, allowing us a complete overview of the viable pMSSM interpretation of the Fermi GeV 
excess. 

The paper is organized as follows: in section 2 we briefly present the well-known frame¬ 
work of the pMSSM and the parameters describing the model. In section 3 we describe the 
experimental set-up of the global fit and the implementation of the joint likelihood. We 
present the results of the parameter scan in section 4. We then discuss the implications for 
future direct and indirect detection experiments in section 5.1 and the prospects for detection 
at the LHC run II in section 5.2. Finally, we summarize our conclusions in section 6. 

2 The theoretical framework: the pMSSM 

We here study the pMSSM [29], in which the number of free parameters can be reduced to 
19, given the present lack of experimental evidence for SUSY and no experimental indication 
that one requires the full freedom of a 22-dimensional pMSSM at present. 

In this model, the lightest supersymmetric particle is the lightest neutralino, %, a com¬ 
bination the neutral electroweak gauginos and higgsinos fields. The neutralino is one of 
the most well-motivated particle DM candidates since it is neutral, stable over cosmological 
timescales and can lead naturally to the correct DM relic abundance in the early Universe. 
In what follow, we assume that the neutralino is the particle DM candidate although we do 
not impose it to fully account for the DM relic abundance as measured by Planck but we 
allow for subdominant contributions to the DM content (cf. section 3). 

We assume first and second generation mass universality, separately in the lepton and 
quark sectors (table 1). The trilinear couplings of the sfermions enter in the off-diagonal parts 
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pMSSM parameters and priors 


Flat 

priors 

Log priors 


Mi [TeV] 

(-5, 5) 

sgn(Mi) log | Mi|/GeV 

(-3.7,3.7) 

Mi [TeV] 

(0.1, 5) 

log M 2 / GeV 

(2,3.7) 

M 3 [TeV] 

(-5, 5) 

sgn(M 3 ) log M 3 /GeV 

(-3.7,3.7) 

m.£ [TeV] 

(0.1,4) 

log mi/ GeV 

(-1,3.6) 

nn-E [TeV] 

(0.1,4) 

log m E / GeV 

(-1,3.6) 

m L3 [TeV] 

(0.1,4) 

log mj j?j / GeV 

(-1,3.6) 

m E3 [TeV] 

(0.1,4) 

log rriE 3 /GeV 

(-1,3.6) 

m Q [TeV] 

(0.1,4) 

log mQ / GeV 

(-1,3.6) 

m.u [TeV] 

(0.1,4) 

log mu /GeV 

(-1,3.6) 

m,£> [TeV] 

(0.1,4) 

log mu/ GeV 

(-1,3.6) 

"iQ 3 I TeV ] 

(0.1,4) 

log mQ 3 / GeV 

(-1,3.6) 

m U3 [TeV] 

(0.1,4) 

logm/ 7 3 /GeV 

(-1,3.6) 

m D3 [TeV] 

(0.1,4) 

log m e 3 / GeV 

(-1,3.6) 

A t [TeV] 

(- 10 , 10 ) 

sgn(At) log At /GeV 

(-4,4) 

[TeV] 

(- 10 , 10 ) 

sgn(A 0 )log|A fe |/GeV 

(-4,4) 

A r [TeV] 

(- 10 , 10 ) 

sgn(A 0 ) log A r /GeV 

(-4,4) 

M [TeV] 

(-5,5) 

sgn(/r) log p /GeV 

(-3.7,3.7) 

m,A [TeV] 

(0.01, 5) 

log 771 , 4 /GeV 

(1,3.7) 

tan/3 

(2,62) 

tan /3 

(2,62) 

M t [GeV] 

173.2 ± 0.87 [33] (Gaussian 

prior) 

po [GeV/cn 

i 3 ] 0.4 ±0.1 [34] (Gaussian prior) 


Table 1 . pMSSM parameters and top mass value used in this paper and the prior range for the 
two prior choices adopted in our scans. “Flat priors” are uniform on the parameter itself (within 
the ranges indicated), while “Log priors” are uniform in the log of the parameter (within the ranges 
indicated). 


of the sfermion mass matrices. Since these entries are proportional to the Yukawa couplings of 
the respective fermions, we approximate the trilinear couplings associated with the first and 
second generation fermions to be zero, while the parameters At, A and A r represent the third 
generation trilinear couplings. In our set-up, the Higgs sector is fully described by the ratio 
of the Higgs vacuum expectation values tan/3, the higgsino mass parameter p, and the mass 
of the pseudoscalar Higgs rri a , which are more directly related to the phenomenology of the 
model. This 19-dimensional realization of the pMSSM encapsulates all phenomenologically 
relevant features of the full model that are of interest for DM and collider experiments. The 
model parameters are displayed in table 1, along with their prior ranges. All of the input 
parameters are defined at the SUSY scale 

3 The experimental setup 

We implement experimental constraints with a joint likelihood function, whose logarithm 
takes the following form: 

In £ joint = In £gce + In Tew + In T B / D ) + In Cn h 2 

x (3.1) 

+ In £lux + In Tic + In Tffiggs + ln ^SUSY, 
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where £qce represents the Fermi GeV excess, Aw electroweak precision observables, A(D) 
B and D physics constraints, A x /i 2 measurements of the cosmological DM relic density, 
Aux (Ac) direct (indirect) DM detection constraints and Aiggs (Ausy) Higgs (sparticles) 
searches at colliders. We discuss each component in turn: 

Ace= For the Fermi GeV excess likelihood we follow the treatment in [15], to account 
astrophysical uncertainties. In addition, we include uncorrelated 10% uncertainties as DM 
modeling systematics, following ref. [32]. We marginalize the likelihood function over the un¬ 
certainties in the Galactic centre J-value, which is assumed to follow a log-normal distribution 
with mean log 10 J/(GeV 2 cm -5 sr -1 ) = 23.29 and variance A log 10 J/(GeV 2 cm' 5 sr -1 ) = 
0.37 [17]. When the predicted neutralino relic density, Q x , is smaller than the Planck mea¬ 
surement, Ddmj we follow [35] and adopt the so-called “scaling Ansatz”, i.e. we assume that 
the local ratio of neutralino ( p x ) to total DM densities (pdm) is equal to that for their cosmic 
abundances: 

Px/P dm = H x /Ddm = f x - (3-2) 

Aw : This implements constraints from Z-pole measurements at LEP [36]. We include the 
constraint on the effective electroweak mixing angle for leptons sin 2 9 e g, the total width of 
the Z-boson Tz, the hadronic pole cross-section (x 9 a(i , as well as the decay width ratios R/, 
R® and the asymmetry parameters Ai, Ab, A c and A°p B , A° B C B . In addition, we also use the 
measurement of the mass of the W boson mw from the LEP experiment [36]. We apply a 
Gaussian likelihood for all of these quantities, with mean and standard deviation as given in 
table II of [37]. 

A(D) : The flavor observables related with B and D physics considered are BR(B —> Xg'y), 

Ram B s , bwb^u)s M ’ A °-> BR ( D s ->• TZ/), Ri23, BR(B s p + p~) and BR(B d p + p~). 
For all we apply a Gaussian likelihood and for most of them we use the measurements 
following table II of [37]. The experimental values assumed for BR{B S —> fi + n~) and 
BR(B d —> are (2.9 ± 0.8) x 1CD 9 and (3.6 ± 1.55) x 1CU 10 (theoretical uncertain¬ 

ties included [38]) [39] respectively. 

A x /i 2: We include the Planck Cosmic Microwave Background data constraint on the DM 
relic abundance as an upper limit, to allow for the possibility that neutralinos are a sub¬ 
dominant DM component. We follow the formalism in the Appendix of [40], using as central 
value the result from Planck temperature and lensing data D x h 2 = 0.1186T0.0031 [41] with a 
(fixed) theoretical uncertainty, t = 0 .012, to account for the numerical uncertainties entering 
in the calculation of the relic density. 

Aux : For DM direct detection we include upper limits from the LUX experiment [42], 
as implemented in the LUXCalc code [43], including both the spin-independent and spin- 
dependent cross-sections in the event rate calculation. We adopt hadronic matrix elements 
determined by lattice QCD [44, 45]. 

Ac: This implements conservative upper limits on the proton spin-dependent cross-section 
from the IceCube detector in its 79-string configuration [46] (IC-79). Comparable - if 
slightly weaker - limits have been set for the WW channel by Super-Kamiokande [47] and 
ANTARES [48]. The most stringent constraint is for the case where WIMPs annihilate ex¬ 
clusively to WW pairs. Since the neutrino spectrum generated by Z bosons is similar to the 
W bosons we apply this constraint whenever the combined branching fraction to WW and 
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ZZ is above 95%. In that case the likelihood is a step function smeared with half a Gaussian 
(as in Eq. (3.5) of [49]) to account for theoretical and experimental uncertainties that we set 
to be 50% of the predicted value. 

£ Higgs' The likelihood for the Higgs searches has two components: the first implements 
bounds obtained from Higgs searches at LEP, Tevatron and LHC via HiggsBounds [50], 
which returns whether a model is excluded or not at the 95% CL. The second component 
constrains the mass and the production times decay rates of the Higgs-like boson discovered 
by the LHC experiments ATLAS [51] and CMS [52], We use HiggsSignals [53] assuming a 
theoretical uncertainty in the lightest Higgs mass calculation of 2 GeV. 

Tsusy: SUSY searches constraints at LEP and Tevatron follow the likelihood used [49]. 
We have imposed the strict constraints from a large number of searches for supersymmetry 
at the LHC. The branching ratios of the sparticles have been calculated with SUSYHIT 
1.5 [54]. We have generated the hadronized event samples with Pythia 8.2 [55] and have 
employed the NNPDF 2.3 parton distribution functions [56]. The generated events are passed 
on to CheckMATE [57], which is based on the fast detector simulation Delphes 3.10 [58]. 
CheckMATE tests if the model point in question is excluded or not at 95% confidence level 
by comparing to current experimental searches at the LHC for supersymmetry in the relevant 
hadronic and leptonic final states with large missing transverse momentum. We assign a 0 
log-likelihood Tsusy = 0 if the point passes all constraints, and exclude it if it fails any of 
them. 

We have only included observables we consider robust in order to be conservative. For 
instance, we have dropped the electroweak precision observables R® and A°p b B in the fit 
because it is unclear whether the large deviations of 2.5<r that are observed with respect to 
the SM predictions are due to unknown systematic uncertainties or to new physics. The 
experimental status of the magnetic anomaly of the muon, = \{g — 2)^ remains unclear in 
the face of persistent discrepancies in the determination of the hadronic vacuum-polarization 
diagram using either e + e - or the hadronic r—decay data. We do not include DM searches 
in dwarf spheroidal galaxies by the Fermi -LAT [59], as constraints are given only for 100% 
branching ratio into final states. 

Finally, we refer to [37] for details about how the SUSY spectrum and observables are 
computed. 

We use the MultiNest [60] algorithm as implemented in SuperBayeS-v2.0, to perform a 
global fit of the pMSSM parameter space, including all the data in eq. (3.1), excepting the 
SUSY searches at the LHC. This is because the LHC searches evaluation is computationally 
too expensive to be performed on-the-fly. Our scans were run using both log and flat priors 
to ensure a complete coverage of the parameter space, gathering ~ 10 6 samples from ~ 10 8 
likelihood evaluations. Samples have been thinned by a factor of 10, focusing our search to 
regions of the parameter space that were not clearly ruled out by LHC-Run I constraints. This 
produced 10 5 representative samples to which the LHC SUSY searches have been applied. 
The ensuing ~ 10 1 samples that pass LHC run I constraints are displayed in figure 1. 


4 Results 

Our global fits identify two distinct viable solutions in the pMSSM parameter space (figure 1): 
the first exhibits a WIMP mass of ~ 80 — 100 GeV, with the neutralino annihilating to WW 
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Figure 1. 2D map of the p -values of our fit, showing the annihilation (top left), the spin-independent 
(top right) and the spin-dependent (bottom panel, left for neutrons scattering and right for protons) 
cross-sections vs neutralino mass. The color-bar represents the p-value from the global fit. The yellow 
overlay highlights points that are within 2er of the Planck relic density. Red crosses indicate the 
best-fit points in the two islands. The Fermi dwarfs limit for the W + W~ channel [59] is plotted for 
reference only, and it has not been applied. To compare with Fermi dwarfs data, the annihilation cross- 
section needs to be rescaled by a factor which would suppress the signal well below current limits. 
The spin-dependent and spin-independent cross-sections have been multiplied by f x = and 

9x = Po/0-3 GeV/cm 3 to facilitate comparison with current and future limits (LUX [42], XenonlT 
and a multi-ton liquid Xe detector with 500 tx yr exposure [61]). In the bottom right panel, we 
display the IC-79 limit [46] used in our analysis. 


with a 95% branching ratio. The second solution has a larger neutralino mass, ~ 180 — 200 
GeV, and 87% it annihilation final states. The overall best-fit point is in the WW region 
has —21n£ = \ 2 = 122.0. This is for a fit with 21 free parameters, and 125 Gaussian data 
points (we do not include limits as their \' 2 is normalized to 0 whenever the constraint is 
satisfied), so we adopt 104 degrees of freedom 1 . Our best fit thus has a p- value of 0.11 versus 

1 We emphasize that the calculation of the number of actual degrees of freedom is not trivial. One would 

have to consider the number of active data points, as well as the number of effective parameters in the model. 
This can only be done properly via extensive Monte Carlo simulations of the data. The simple counting 
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Figure 2. Contribution to the overall x' 2 for the two best-fit points, grouped by type of observable (see 
section 3 for details). The pulls have been normalized by the number of data points in each group, 

N, to facilitate a visual comparison. The dashed bars show the Galactic Centre Fermi likelihood 
contribution when the 10% theoretical uncertainty is neglected, which would degrade the p-values to 

O. 023 (for the m x = 88.3 GeV solution) and 0.008 (for the m x = 188 GeV model). 

a x 2 = 127.6 and a p-value of 0.06 for the it solution. 

It is important to notice that including theoretical uncertainties in the GC fit is crucial 
in achieving reasonable p-values. The 10% theoretical uncertainty advocated in [32] is a 
reasonable reflection of the differences in the predicted spectra between current numerical 
codes. However, in absence of such an uncertainty, the quality of our global fit would degrade 
to 0.023 and 0.008 for the two best-fit points, respectively. 

The contribution to the overall x 2 for the two best-fit points from different types of 
observables are plotted in figure 2. The pulls have been normalized by the number of data 
points in each group, IV, to facilitate a visual comparison. We notice that the x 2 P er data 
point is distributed fairly evenly across observables. There is a slight preponderance in the x 2 
contribution coming from the GC Fermi fit (with N = 24 bins), which is exacerbated if one 
neglects the 10% theoretical uncertainty in the DM spectra (dashed bars). The contribution 
to the pull from the LUX likelihood comes almost exclusively from the SD neutron cross- 
section limit, as the SI constraint is easily satisfied by our best-fit points. 

Let us now analyze in more detail both type of solutions. 

In the WW region, model points providing a better fit have a neutralino mostly bino- 
like (~ 80 - 90%) with a similar fraction of both wino and higgsino. Besides we find points 
in which neutralinos can be dominantly higgsinos with a bino fraction as small as ~ 10%. 
Those provide a worse fit though, because a large higgsino composition, basically, implies a 
large annihilation cross section which drops the relic density below the Planck limit leading 
to a tension with the Fermi GeV excess due to the scaling Ansatz given in eq. (3.2). 

argument we adopt is meant to be representative of what one would get in the simplest scenario. 
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Analogously to Ref. [32], also in the present work, we find a WW solution with a bino- 
higgsino neutralino composition of about 84 - 92 GeV. However, in our work, this solution 
has a slightly worse x 2 than in Ref. [32] . The reason for this is that we adopt slightly different 
values of the form-factors for the computation of the SD cross-section. Since such a WW 
solution is now just below the WW IceCube limits, it is punished a bit in the likelihood. 

Typically, since third generation squarks enter into the SI cross section at loop level 
through LSP-gluon effective interactions [62], their contribution can be comparable to tree- 
level effective interactions mediated by squarks of the first two families when they are light. 
In this type of solutions, because current constraints allow sbottoms of a few hundred GeV, 
their contribution can be sizable and, indeed, cancels out the Higgs exchange contributions. 
This effect allows to relax the tension with LUX data when Higgs exchange contributions are 
large. 

In terms of the impact of LHC Run-I data, we notice that first and second generation 
squarks as well as gluinos are decoupled. The stops and the heavier sbottom mass eigenstates 
are also kinematically inaccessible at the LHC Run-I energies. However, the lighter sbottom 
eigenstate with a mass around 400 GeV can be produced at the LHC at a considerable rate 
but the sbottoms evade detection from third generation searches due to complicated cascade 
decays. The sparticles in the electroweak sector are relatively light. There is a large mass 
splitting between the SU(2) doublet and singlet sleptons. The SU(2) doublet sleptons are 
light with masses around 250 GeV and narrowly evade detection in searches for two lepton 
and large missing transverse momentum final states. The lightest wino-like chargino and the 
second lightest neutralino escape detection since they are almost mass degenerate with the 
bino-like neutralino. The production rate of the higgsino eigenstates is too small to yield an 
observable signal at the LHC run I. 

Results about the annihilation cross-section are shown in the top left panel of figure 1. 
One can see that the points with a better fit exhibit (av) ~ 1CU 26 cm 3 /s, consistent with 
the results found in [15]. We also show the Fermi dwarfs limit for the W + W~ channel [59], 
but we emphasize that this limit has not been applied in the fit. In order to compare with 
the constraint coming from Fermi dwarfs observations, the annihilation cross-section needs 
to be rescaled by a factor / 2 (to account for the possibility of sub-dominant dark matter relic 
density, which translates according to our Ansatz in a correspondingly reduced local density 
in the dwarfs), which would suppress the signal well below current and future limits from 
dwarfs. 

Regarding DM direct detection, the top right panel of figure 1 shows the spin-independent 
(SI) cross-section versus neutralino mass plane. In order to facilitate the visual comparison 
of our pMSSM models with existing and future limits, we have rescaled the theoretical cross- 
section by a factor f x (to account for models where the neutralino does not make up all of 
the cosmological relic density) and a factor g x = po/0.3 GeV/cm 3 . This accounts for the fact 
that the local density we have used to predict the number of counts for LUX, po, is a nuisance 
parameter which is generally different from the value assumed by the LUX collaboration in 
deriving their limit [42], namely 0.3 GeV/cm 3 . 

The points that appear above the nominal 95% exclusion limit from LUX cannot be 
excluded because of a combination of effects: (i) our LUX likelihood is slightly less stringent 
than what has been published by the LUX collaboration (and depicted in figure 1), (ii) our 
likelihood function for LUX allows for values of the cross section above the 95% limit to be 
included (albeit penalized by a smaller likelihood value), and (iii) the global likelihood func- 
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tion can - to an extent - compensate for a poor fit to LUX data by gaining an improvement 
from other data sets. 

Since squarks are heavy, the contribution coming from the exchange of a CP-even Higgs 
dominates and therefore the SI cross-section scales as oc |A r nA r 13 / 14 | 2 . Because, as noticed 
above, the higgsino fraction in this region is not negligible, it can be large. Indeed, there are 
model points above the LUX limit allowed by the scaling Ansatz we apply to the local dark 
matter density and to some extend due to the fact that we vary the local dark matter density. 
Another interesting feature is that the SI cross-section spans down to 10 14 pb. That is 
possible because the heavy CP-even Higgs contribution can be sizable and cancellations with 
the lightest Higgs channel might occur [63]. 

In the bottom panels of figure 1 we display the spin dependent (SD) cross-section for 
scattering off neutrons/protons (left/right panels). While in the case of SI interactions the 
contributions for proton and neutron are comparable, the SD cross-sections may differ sig¬ 
nificantly. However, we find a tight correlation between the SD cross-sections for scattering 
off neutrons and those off protons in our model points. The SD cross-section is dominated 
by the exchange of a Z boson and therefore the SD cross-section is largely determined by the 
higgsino content of the neutralino and likewise for the SI cross-section it can be sizable in 
this region as it can be seen in both panels. 

In terms of dark matter direct detection experiments, at present, LUX data represents 
the strongest constraint on the SD-neutron scattering cross-section because the Xenon con¬ 
tains neutron-odd isotopes therefore we overlay the LUX constraint properly rescaled as for 
the SI case. For the SD-proton scattering cross-section, IC-79 represents the strongest cur¬ 
rent constraint for the particular case when the neutralinos annihilate to a W + W _ final state 
so we show the IC-79 90% CL upper limit [46] , and we have rescaled the value of the SD 
cross-section by a factor f x . This assumes that equilibrium between capture and annihilation 
is reached in the Sun, which is a good assumption for the bulk of the models shown here. 
In fact, in the region where the neutralino annihilates mainly to WW and to ZZ to a lesser 
extend, the IC-79 limits apply here and disfavors a large number of model points. 

One can see that there are model points above LUX and IC exclusion lines. In those, the 
higgsino component of the neutralino is dominant over the gaugino one leading to a large SD 
cross-section for both neutrons and protons. Those points still provide a reasonable fit to the 
data due to two effects, first because neutralinos with a large higgsino component yields to 
a relic density sensible below the Planck bound and therefore the scaling Ansatz applies and 
second because the local dark matter density is a nuisance parameter in our analysis. Beside 
there is a small fraction of points above the nominal IC-79 WW limit with a branching ratio 
to WW and ZZ low enough to evade the stringent IC-79 bound. 

In the it region, points providing a better fit have a neutralino dominantly bino-like 
(~ 90%) with a ~ 10% of higgsino. Those points have the characteristic that the neutralinos 
annihilate to top quark pairs via an exchange of a right-handed stop which is relatively heavy 
(~ 1 TeV). This is possible because on top of the non-helicity suppression the neutralino- 
stop-top coupling component, which is proportional to the top quark Yukawa coupling, is 
sizeable due to the non-vanishing higgssino fraction of the neutralino. We also found points 
in which the right-handed stops are light (~ 300 GeV) being almost bino-like (also found 
in [32]). However those provide a worse fit to the BR(B —> Xg'y) data because higgssino- 
stop loops have a sizable positive contribution which leads to values above the experimental 
constraint. 
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Annihilation into it through a t-channel stop exchange requires, for bino-like neutralinos, 
stops with masses of a few hundred GeV. On the other hand, relatively light stops are not 
able to lift the tree-level Higgs mass to fit a 125 GeV Higgs. In order to enhance the DM 
annihilation cross sections to match the Planck measurement, the neutralino coannihilates 
with sneutrinos which have to be about the same mass as the DM particle (~ 200 GeV). 
This induces the splitting in the left/right sleptons spectrum. 

These benchmark points are characterized by left handed slepton next-to-lightest super- 
symmetric particle with masses above 200 GeV. The higgsinos have masses around 260 GeV. 
Since the resulting mass differences between the bino-like LSP and the sleptons/higgsinos 
are small, the final states are rather soft and thus the detection is suppressed in events with 
dilepton or trilepton final states and large missing transverse energy momentum. The pro¬ 
duction rate of the sbottom is quite suppressed and hence avoids detection. The spectrum 
of the remaining supersymmetric particles is decoupled. 

The phenomenology of the model points in the it region in terms of DM detection is 
similar to the WW one. The main difference is in that the IC-79 limits does not apply 
and therefore larger SD cross-sections are possible. The main difference is that the higgsino 
composition is not as large as in the WW type of solutions and therefore the SD cross sections 
for neutron interactions are below the LUX current sensitivity. 

Lastly, in figure 3, we show the spectrum of the Fermi GeV excess together with the 
systematic uncertainties associated with the galactic diffuse emission modeling [15]. We 
compare the data with the spectra of the pMSSM model points giving the best global p- 
value in the two regions identified in Figure 1. 

It is apparent from Figure 3 that the best-fit DM spectra are systematically offset from 
the mean values of the Galactic center excess spectrum (gray dots and boxes), by about 1 
to 2 sigma, and do not provide a good fit to the data at first sight. However, since the 
systematic astrophysical uncertainties, indicated by the gray boxes (±lcr) are correlated , this 
still provides an acceptable fit to the data. To illustrate this point, we show with black dots 
and error bars the excess spectrum where we moved all data points systematically down, 
according to the freedom allowed by the covariance (the error bars show now statistical 
errors only). Together with the 10% uncorrelated systematic modeling uncertainty that we 
adopted for the DM signal, this provides a reasonable fit to the data, with p-values, whereas 
without the DM signal modeling uncertainties, the p -values would be prohibitively small (see 
Figure 2 above). 

5 Discovery potential: dark matter detection experiments and LHC run 

II 

5.1 Implications for direct and indirect dark matter searches 

Generally, spin-dependent and spin-independent scattering cross sections are driven by the 
higgsino content of the neutralino. Therefore, as explained above, the sizable higgsino frac¬ 
tion of points in the two best-fit regions imply large SI cross-sections that makes the direct 
detection prospects promising, although the SI cross-section range spans down to ~ 10“ 14 
pb due to cancellations with the lightest Higgs. In figure 1, top right panel, we also display 
the projected sensitivity limit (defined as the 90% CL exclusion limit) for the XenonlT ex¬ 
periment and an hypothetical liquid Xe detector with 500 tx yr exposure [61]. The latter 
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Figure 3. Spectral energy distribution of the Fermi -GeV excess data. Grey dots show Fermi GeV 
excess mean values w.r.to background model variations, together with associated systematic uncer¬ 
tainties {grey boxes). Black dots represent the excess for a variation of the Galactic diffuse emission 
contribution (within its systematic uncertainty). The solid lines show the prediction for the pMSSM 
models that give the best fits in the two regions of figure 1 and have dominant annihilation channel 
into W bosons or top quarks {magenta lines in the left and right panel, respectively). The green lines 
indicate the adopted 10% systematic uncertainties in the spectra. Furthermore, A denotes the re¬ 
quired boost-factor with respect to a generalized Navarro-Frenk-White DM density profile with inner 
slope 7=1.26. The p denotes the p -value of the global fit (including all data). 


experiment essentially saturates the ultimate detection floor set by coherent neutrino scatter¬ 
ing [64]. XenonlT data will be crucial in discovering or firmly ruling out models belonging to 
the tt island and will probe a significant fraction of the parameter space preferred by the first 
type of solutions. In figure 1, bottom left panel, we overlay the projected 90% exclusion limit 
for XenonlT [61]. We find that XenonlT will be able to prove the entirety of the SD neutron 
scattering cross-section parameter space favoured by our models. 2 A multi-ton experiment 
with 500 tx yr exposure would reach sensitivities a couple of order of magnitudes smaller 
than the smallest SD neutron cross-sections found in our scan. 

The bottom right panel of figure 1 shows that our best-fit points easily evade the con¬ 
straint set by IC-79 on the SD proton cross-section in the implementation we adopted in 
this paper. However, an event-level implementation of the likelihood (including the events’ 
energies [66]) would increase the constraining power of the IC-79 limit, to the point that 
some of the surviving models could be probed [67]. 

Finally, as for indirect detection, the preferred parameter space is mostly out of reach 
even for the future 10-yr Fermi analysis of dwarf spheroidal galaxies [68], where an further 
improvement of current sensitivities by a factor 2-3 can be expected. Note that points 
with a large annihilation cross section as shown in the top right panel of figure 1 usually 
correspond to suppressed relic densities, making these points hard to detect due to the / 2 
factor suppression of their signal. 
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Figure 4. 2D map of the p-values of our fit, showing the lightest squark mass of the first and 
second generation (top left), the lightest stop (top right), the lightest sbottom (middle left), the 
gluino (middle right), the lightest srnuon (bottom left) and the lightest chargino (bottom right) vs 
the neutralino mass. The yellow overlay shows points within 2<r of the Planck relic density value. 
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5.2 Prospects for detection at the LHC run II 

The models yielding the best p-values show several interesting properties relevant for LHC 
searches. In the following we briefly discuss the discovery prospects depending on the pro¬ 
duced SUSY particles. Sparticle mass distributions and MSSM parameters relevant for LHC 
searches are shown in figure 4. 

Light squarks : In the top left panel of Fig. 4 we plot the lightest first/second generation 
squark mass vs the neutralino mass. The first and second generation squarks have masses 
> 1400 GeV, i.e. above the usual run I constraints. About 70% of the models have squark 
masses below 2000 GeV ’. The upcoming run II searches for light squarks will exclude many 
scenarios. In almost all models the left- and right-handed squarks have different branch¬ 
ing ratios and tend to decay to the heavy neutralino and chargino states. Cascade decays 
including W, Z and Higgs bosons are common. 

Stop: In the top right panel of Fig. 4 we plot the stop mass vs the neutralino mass. Some 
models have light stops with masses down to 200 - 300 GeV decaying to chargino and a b-jet. 
The neutralino has a mass of around 95 GeV. These models are not excluded by current LHC 
searches [69]. Another interesting region also found in [32] has a stop mass of around 200 - 
220 GeV and a mass of the lightest neutralino around 180 GeV. A slight excess in the ATLAS 
data prevents exclusion with run-1 data in this region [69]. These solutions will likely be 
tested with early run II analyses. Other solutions yield much heavier stop masses decaying 
predominantly into the heavier neutralino and chargino states. Dedicated searches for such 
decays are important. 

Sbottom: In the central left panel of Fig. 4 we plot the sbottom mass vs the neutralino mass. 
Several model points have a sbottom mass as low as 400 GeV. The points are not excluded 
in our procedure due to multi-step cascade decays involving heavy neutralinos. 

Typically, the lighter sbottom state typically has masses around 400 GeV while the 
lightest neutralino is mostly bino-like with a mass around 90 GeV. The second lightest wino- 
like neutralino lies around 107 GeV. The two heavier neutralinos are higgsino dominated 
states and have masses around 250 GeV. 

If the sbottoms predominantly decayed into a bottom quark and the lightest neutralino, 
these benchmark points would clearly be excluded. However, the bottom squark decays into 
all four neutralino as well as the lighter chargino mass eigenstates with comparable rates. The 
decay modes of the neutralino and chargino mass eigenstates are relatively complex and hence 
the limits from simplified sbottom searches do not apply. For instance, the second lightest 
neutralino eigenstate has large hadronic three body decay modes into the lightest neutralino 
via off shell Z bosons. Moreover, the second lightest neutralino can radiatively decay into 
a photon and the lightest neutralino. The corresponding lighter chargino eigentstate has 
relatively large leptonic three body decay modes via off shell W bosons. Finally, the two 
heaviest higgsino dominated neutralino mass eigenstates mainly decay into electroweak gauge 
bosons and the lighter electroweakino states. As a consequence, many events have leptons in 

2 We note that the XenonlT exclusion limits in [61] are obtained by applying a scaling factor derived from 
the comparison between SI and SD results of XenonlOO and result to be slightly stronger than the limits 
quoted in ref. [65], which instead adopts a 60 times improvement w.r.to XENON100 (still XenonlT will be 
able to entirely probe our identified parameter space). 

3 The expected reach of the ATLAS and CMS experiments for full MSSM models is not estimated. We 
compare with mass scales which are likely excluded with HL-LHC data for many MSSM scenarios. 
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the final state which are vetoed in the searches for direct sbottom production. In addition, 
the higher final state multiplicity tends to soften the net missing transverse momentum 
distribution compared to scenarios with direct sbottom decays into a bottom quark and the 
LSP. As a result fewer events pass the selection cuts of the relevant sbottom searches. We 
explicitly tested those light sbottom scenarios against experimental searches at the LHC with 
the computer tool CheckMATE and confirmed that those points were allowed. 

As seen also for the light squarks about 70% of the models would be excluded with a 
sbottom limit of « 1000 GeV. 

Gluinos: In the central right panel of Fig. 4 we plot the sbottom mass vs the neutralino 
mass. Gluinos have masses > 1600 GeV. An upcoming early run II exclusion on gluinos with 
a mass up to 2500 GeV would exclude about 15% of these models. Run II searches for first 
and second generation squarks will thus likely be more constraining. 

Sleptons: In the bottom left panel of Fig. 4 we plot the lightest smuon mass vs the neutralino 
mass. The lightest smuons found in the best fit models have masses < 400 GeV in about 
60% of the best fit points. This makes searches for smuons in run II very sensitive to these 
solutions. 

Chargino/neutralino: In the bottom right panel of Fig. 4 we plot the chargino mass vs 
the neutralino mass. Due to the GeV excess likelihood several neutralinos and charginos are 
typically light, the lightest having a mass fixed to 80 — 100 GeV for the WW solutions and 
180 — 200 GeV for the it solutions. The higgsino component in the it and WW solutions 
typically involves that n is only slightly larger than these mass scales, leading to 2 more 
neutralinos and the light chargino at masses around 100 GeV {WW) or 200 - 300 GeV {it). 
These states are often mass compressed with the lightest neutralino which makes the solutions 
evade LHC chargino/neutralino searches so far. Dedicated chargino/neutralino searches will 
have sensitivity to most models, e.g. by a mono-jet and soft lepton search as proposed also 
in ref. [32], The Wino mass scale is quite unconstrained and lies between 100 GeV and 1.5 
TeV. The Wino will decay to lighter states yielding final states with Z, W or Higgs bosons. 

Heavy Higgs: About 50% of the best models have m\ < 800 GeV making searches for 
heavy Higgs bosons very sensitive. Several chargino and neutralino states are light and have 
a large coupling to A/H/H . Consequently heavy Higgs decays to charginos and neutralinos 
can have huge branching ratios up to 30% competing with top and bottom decays. Dedicated 
searches for heavy Higgs bosons decaying into final states with WjZ/h with missing transverse 
momentum would help to constrain these scenarios. 

6 Conclusions 

In this paper we addressed the issue of finding model points in the pMSSM that might 
explain simultaneously the large set of independent data we gathered from astrophysics, 
cosmology and high-energy particle physics. We showed that no tension exists between 
currently available particle physics constraints and the interpretation of the Fermi GeV excess 
in terms of dark matter annihilation in the framework of the pMSSM. Furthermore, we found 
evidence for two regions that are able to explain the gamma-ray data, while being consistent 
with other various experimental constraints: (a) a first region where the neutralino is mostly 
bino-like and the dominant annihilation channel today is 95% into W bosons pairs and (b) a 
second region where the the annihilation into top-quark pairs dominates and the neutralino 


- 14 - 


is again mainly bino-like. We showed that these models are very appealing since they will be 
soon in the reach of the next generation of direct detection experiments - XenonlT will probe 
the entirety of the best-fit regions thanks to its sensitivity to spin-dependent neutron cross- 
section - and of the LHC run II, in particular through searches for charginos and neutralinos, 
squarks and light smuons. 
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